data <- read.table(file.choose(), header=T, sep="\t")
attach(data)

outFile <- file.choose(new=T)


CalcCoeficient <- function(Ave, Dis, chr){
	 Ave <- Ave[Type==chr]
	 Dis <- Dis[Type==chr]
	
	 FITTING_MAX <- max(Dis[Dis <= (max(Dis) - 1000000)])
	 if(FITTING_MAX < 2.3e6){
		FITTING_MAX <- 2.3e6
	 }
	 y <- Ave[Dis < FITTING_MAX & Dis > 0 & Ave != 0]
	 x <- Dis[Dis < FITTING_MAX & Dis > 0 & Ave != 0]
	 a0 <- min(y)
	 c0 <- log((y[1]-a0)/(y[10]-a0))/(x[10]-x[1])
	 b0 <- (y[10]-a0)/exp(-x[10]*c0)/11
	 e0 <- log((y[1]-a0-b0*exp(-c0*x[1])) / (y[3]-a0-b0*exp(-c0*x[3]))) / (x[3] - x[1])
	 d0 <- (y[1] - a0 - b0 * exp(-c0*x[1])) / exp(-e0*x[1])
	 browser()
	 #fit <- nls(y ~ a + b*exp(-c*x) + d*exp(-e*x), start=c(a=a0, b=b0, c=c0, d=d0, e=e0))
	 fit <- nls(y ~ a + b*exp(-c*x), start=c(a=a0, b=b0, c=c0))
	 a <- fit$m$getPars()["a"]
	 b <- fit$m$getPars()["b"]
	 c <- fit$m$getPars()["c"]
	 d <- fit$m$getPars()["d"]
	 e <- fit$m$getPars()["e"]

	 ymax <- max(Ave)
	 plot(Dis, Ave, ylim=c(0,ymax), xlab='', ylab='average frequency')
	 par(new=T)
	 new <- data.frame(x=Dis)
	 prediction <- predict(fit, new)
	 plot(Dis, prediction, ylim=c(0,ymax), type='l', col='blue', axes=F, xlab='', ylab='')
	 legend("topright",lty=1, lwd=2, col="blue", paste(sprintf("%.3f",a0), "+",
		sprintf("%.2f",b0), "*exp(-", sprintf("%.2e",c0), "*x)+",
		sprintf("%.2f",d0), "*exp(-", sprintf("%.2e",e0), "*x)", sep=""))
	 t1 <- paste(chr, "a", a0, sep="\t")
	 t2 <- paste(chr, "b", b0, sep="\t")
	 t3 <- paste(chr, "c", c0, sep="\t")
	 t4 <- paste(chr, "d", d0, sep="\t")
	 t5 <- paste(chr, "e", e0, sep="\t")
	 cat(t1, t2, t3, t4, t5, sep="\n", append=TRUE, file=outFile)
}

par(mfrow=c(3,1), oma=c(1,1,4,1), mar=c(3,4,0,0))
c1 <- CalcCoeficient(Average, Distance, "chr1")
c2 <- CalcCoeficient(Average, Distance, "chr2")
c3 <- CalcCoeficient(Average, Distance, "chr3")

t6 <- paste("inter", "chr1_chr2", Average[Type=="chr1_chr2"], sep="\t")
t7 <- paste("inter", "chr1_chr3", Average[Type=="chr1_chr3"], sep="\t")
t8 <- paste("inter", "chr2_chr3", Average[Type=="chr2_chr3"], sep="\t")
cat(t6,t7,t8, sep="\n", append=TRUE, file=outFile)


